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We study the kinetics of the nematic-isotropic transition in a two-dimensional liquid crystal by 
^— «j ' using a lattice Boltzmann scheme that couples the tensor order parameter and the flow consistently. 

I Unlike in previous studies, we find the time dependences of the correlation function, energy density, 

, ^ ■ and the number of topological defects obey dynamic scaling laws with growth exponents that, within 

^ ' the numerical uncertainties, agree with the value 1/2 expected from simple dimensional analysis. 

We find that these values are not altered by the hydrodynamic flow. In addition, by examining 
shallow quenches, we find that the presence of orientational disorder can inhibit amplitude ordering. 
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I. INTRODUCTION 



■ The phase orderingkinetics of liquid crystal systems undergoing nematic-isotropic transitions has attracted consid- 
erable experimental [0-^ and theoretical [§-0 interest. One of the reasons is that nematic liquid crystals provide an 
experimentally accessible system with continuous symmetry that, unlike systems with discrete symmetry, can display 
stable topological defects. During phase ordering these defects interact and annihilate and it is believed that many 
of the universal properties of the late stage kinetic growth can be explained in terms of the defect dynamics. This 

^ property is also shared by the well studied 0{N) model. However here the continuous symmetry belongs to a different 
homotopy class (the 0{N) model lacks the inversion symmetry present in the director field of nematic liquid crystals). 
It is not yet completely clear if such a difference in the symmetry can modify the ordering dynamics of the nematic 
liquid crystal. 

The aim of this paper is to study the kinetics of phase separation for two-dimensional liquid crystals under a quench 
iy~j from an isotropic to a nematic phase. Considerable controversy exists as to whether the ordering violates dynamical 
scaling The simplest scaling analysis based on the assumption of diffusive dynamics for the order parameter 
' suggests that any length scale in the system should grow with a power law in time t^/^ However, an examination 
[ of the configuration of the order parameter in experiment (or simulation) clearly shows that the late stage ordering 
1^ ' proceeds by defects moving to annihilate. Simple defect arguments which rely on an assumption of a finite constant 
' friction coefficient for the movement of defects also give t^/^ [Q. A problem arose due to early calculations [§-0 which 
"j^ showed that the friction coefficient diverges logarithmically with system size. This brought into question whether the 
t^/^ behavior from the simple scaling analysis would be found. Simulations of the XY model [|| and later of tensor 
i' , models of liquid crystals |^J^ reinforced this view when they failed to measure the i^^^ behavior and, in fact found 
'■^ ' exponents which appeared to be decreasing away from 1/2 at late times However, experiments by Pargellis and 
^ [ coworkers ||^ found behavior consistent with the t^^^ power law. In addition, their simulation of the XY model with 
O ■ very large amplitude noise (considerably larger than that found experimentally) agreed with the <^/^ power law. A 
later more rigorous calculation of the friction coefficient for defect dynamics shows that it does not diverge with 
system size, but only with vanishing defect velocity |l^. This suggests that it should be possible to measure the t^/^ 
behavior in a simulation. 

Another important issue which remains essentially unexplored is the extent to which the presence of hydrodynamic 

■ modes coupled to the nematic order parameter can affect the late stage kinetic growth. It is well known that 
hydrodynamic interactions can crucially infiuence phase transition kinetics, for example in simple binary fluid mixtures 
and in gas- liquid systems []l2| , ^ . In liquid crystals the coupling between the local molecular orientation and the 
velocity field influences the dynamics of the liquid crystal in a complicated manner. For example, shear flow will cause 
the molecules to reorient and conversely a reorientation of the liquid crystal may induce a velocity field (backflow). 
To have a complete picture of the kinetic growth of a phase separating nematic liquid crystal it is crucial to include 
the effect of these couplings. 

Lattice Boltzmann approaches have proved very successful algorithms for investigating phase ordering in binary 
fluids. Here we use an extension of the method which models liquid crystal hydrodynamics [ p^JT5| . The equations 
of motion describing liquid crystal hydrodynamics are complex. There are several derivations broadly in agreement, 
but differing in the detailed form of some terms. Here we follow the approach of Beris and Edwards who write the 
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equations of motion in terms of a tensor order parameter Q which describes the orientational distribution function of 
the molecules . This formalism is appropriate here because the motion of defects is explicitly included. Moreover 
both the isotropic and nematic phases can be modelled using the same formalism which is necessary if the dynamics 
of the transition between them is to be followed. 

The paper is organized as follows. In Section II wc summarize the hydrodynamic equations of motion for liquid 
crystals. The lattice Boltzmann scheme used to model these equations is described in Section III. In Section IV we 
introduce the correlation function, the correlation length, the energy density, and the defect separation as the relevant 
quantities that characterize the late time behavior of the phase separating system. A review and discussion of previous 
results is also given. The dynamic scaling behavior of the model is investigated for systems of linear size L = 256, 
and the space of the relevant parameters of the model is explored to choose the best values for more computationally 
intensive simulations. In Section V a careful estimate of the dynamical exponents is performed for systems of linear 
size L — 512. Numerical results for shallow quenches are presented in Section VI. In Section VII we make some 
concluding remarks. 



II. THE HYDRODYNAMIC EQUATIONS OF MOTION 



There are two major differences between the hydrodynamics of simple liquids and that of liquid crystals. First, 
the rod-like shape of the molecules means that they are rotated by gradients in the velocity. Second, the equilibrium 
free energy is more complex than for a simple fluid and this in turn increases the complexity of the stress tensor 
in the Navier-Stokes equation for the evolution of the fluid momentum. In a general formulation of liquid crystal 
hydrodynamics , the continuum equations of motion are written in terms of a tensor order parameter Q which is 

traceless and symmetric and is related to the direction of individual molecules rh by Qap — {'rha'n^ii — ^Sa/s) where the 
angular brackets denote a coarse-grained average. (We shall use Greek indices to represent Cartesian directions and 
assume the usual sum over repeated indices.) The advantage of this approach is that it includes both the isotropic 
(Q = 0) and the nematic (Q ^ 0) phases and allows an order parameter of variable magnitude within the latter. Hence 
it is possible to explore the effect of flow on the phase transition between the two states. Moreover the hydrodynamics 
of topological defects (point defects in two dimensions) is naturally included in the equations. We will study a two- 
dimensional system with the flow confined to the xy-plane. However we will allow the director field to point in any 
direction {x, y and z). As such the tensor order parameter is always a 3x3 matrix. 

The equilibrium properties of the liquid crystal are described by the Landau-de Gennes free energy functional |l^] 

^ = I d'rS^^il - ^)Ql0 - ^Qaf^Qp.Q.a + \{Qlpf + \(dMpxf\^ ■ (1) 

This free energy describes a first-order transition from the isotropic to the nematic phase. In general three elastic 
constants are needed to fully characterize the nematic phase ||l^ but we restrict ourselves to a single elastic constant 
K. This can be shown to be equivalent to having three equal Frank elastic constants when the order parameter is 
uniaxial This simplification is not expected to affect the scaling behavior. 

The order parameter Q is not conserved. It evolves according to a convection-diffusion equation p6| , pj -21 



(St+u- V)Q-S(W,Q) =fH (2) 
where u is the bulk fluid velocity and f is a collective rotational diffusion constant. A generalized form of f is 

^ " (l-|TrQ2)2 

where the Q-dependence in the denominator enhances reorientation for well-ordered systems pO| . Note that in 
previous studies of the kinetics of phase separation for liquid crystals the Q-dependence in the diffusion parameter 
has been always neglected and F = F assumed. 

The first term on the left-hand side of equation (|^) is the material derivative describing the usual time-dependence 
of a quantity advected by a fluid with velocity u. This is generalized by a second term 

S(W, Q) = (D + fl)(Q + 1/3) + (Q + I/3)(D - O) - 2(Q + I/3)Tr(QW) (4) 

where D = (W -I- W^)/2 and O = (W — W-^)/2 are the symmetric part and the anti-symmetric part respectively of 
the velocity gradient tensor W^ji = dpUa- S(W, Q) appears in the equations of motion because the order parameter 
distribution can be both rotated and stretched by flow gradients. 
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The term on the right-hand side of equation (|^) describes the relaxation of the order parameter towards the 
minimum of the free energy in a way analogous to Model A [^2| . The molecular field H which provides the driving 
motion is related to the derivative of the free energy by 

The flow of the liquid crystal of density p obeys the continuity 

dtp + dapua = (6) 

and the Navier-Stokes equation 

pTf 

pdtUa + pUpdpUa = dpTap + BpGap + — 9/3((l - id pPq)5 apd-^U^ + d^Up + dpUa) (7) 

where t/ is related to the viscosity and Pq is the pressure, 

Po = pr-|(VQ)2. (8) 

The details of the stress tensor reflect the additional complications of liquid crystal hydrodynamics with respect to 
simple fluids. There is a symmetric contribution 

Cfap = — Po^Q/3 — iHap ~ dpQ^y (9) 

and an antisymmetric contribution 

Tal3 = QajHjp — HajQjp. (10) 

For the symmetric contribution we are using the form derived by Doi [Q. This is only quantitatively correct in the 
vicinity of the transition with the general form being slightly more complex . We do not expect any qualitative 
differences to result from this difference in the regime in which we operate in this paper. 

III. A LATTICE BOLTZMANN ALGORITHM FOR LIQUID CRYSTAL HYDRODYNAMICS 

We now define a lattice Boltzmann algorithm which solves the hydrodynamic equations of motion of a liquid crystal 
@i ® J and (^. This section may safely be omitted by readers interested in the physical results but not in the details 
of the simulations. 

Lattice Boltzmann algorithms are defined in terms of a set of continuous variables, usefully termed partial distri- 
bution functions, which move on a lattice in discrete space and time p3|] . The simplest lattice Boltzmann algorithm, 
which describes the Navier-Stokes equations of a simple fluid, is defined in terms of a single set of partial distribution 
functions which sum on each site to give the density. For liquid crystal hydrodynamics this must be supplemented by 
a second set, which are tensor variables, and which are related to the tensor order parameter Q. 

We define two distribution functions, the scalars fi{x) and the symmetric traceless tensors Gi{x) on each lattice 
site X. Each fi, Gi is associated with a lattice vector e^. We choose a nine- velocity model on a square lattice with 
velocity vectors = (±1, 0), (0, ±1), (±1, ±1), (0, 0). Physical variables are defined as moments of the distribution 
function 

P = ^fi^ PUa=^fieia, Q = ^Gj. (11) 

i i i 

The distribution functions evolve in a time step At according to 

+ e,At, t + At)- f,{x, t)^^ [Cf^{x, t, {/,}) + Cf,{x + e,At, t + At, {/*})] (12) 

Gi(x + e,At, t + At)- Gdx, t) = ^ [Ccrix, t, {Gi}) + Cg^{x + e,At, t + At, {G*})] (13) 
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The left-hand side of these equations represents free streaming with velocity e^, while the right-hand side is a 
collision step which allows the distribution to relax towards equilibrium. /* and G* are first order approximations to 
fi{x + eiAt,t + At) and Gi{x + eiAt,t + At) respectively. They are obtained from equations ( p^ and (1^) but with /* 
and G* set to fi and G^. Discretizing in this way, which is similar to a predictor-corrector scheme, has the advantages 
that lattice viscosity terms are eliminated to second order and that the stability of the scheme is improved. 

The collision operators are taken to have the form of a single relaxation time Boltzmann equation, together with a 
forcing term 

Cf,{x, t, {f,}) = --{Mx, t) - ft\x, t, {/,})) + MS, t, {/,}), (14) 



CG.(f , t, {G,}) = (G,(f, t) - G^«(f, t, {G,})) + M,(f , t, {GJ). (15) 

tg 

The form of the equations of motion and thermodynamic equilibrium follow from the choice of the moments of the 
equilibrium distributions f^"^ and G^' and the driving terms pi and M^. f^"^ is constrained by 

i i i 

where the zeroth and first moments are chosen to impose conservation of mass and momentum. The second moment 
of f^'' controls the symmetric part of the stress tensor, whereas the moments of pi 

= 0, ^PtBia = djBTafi, '^P^eiae^p = (17) 
i i i 

impose the antisymmetric part of the stress tensor. For the equilibrium of the order parameter distribution we choose 

Gi"^ = Q, ^ Gl'^eia = QUa, ^ G'r'^e^aei|3 = QUaUp. (18) 

i i i 

This ensures that the order parameter is convected with the flow. Finally the evolution of the order parameter is 
most conveniently modeled by choosing 

^M, =fH(Q) + S(W,Q), ^M,e,„ = (^M,)u„. (19) 

i i i 

which ensures that the fluid minimizes its free energy at equilibrium. 

Conditions (|l|)-(^ can be satisfled as is usual in lattice Boltzmann schemes by writing the equilibrium distri- 
bution functions and forcing terms as polynomial expansions in the velocity |p3|| . Taking the continuum limit of 
equations (^2|) and ( [l3| ) and performing a Chapman-Enskog expansion leads to the equations of motion of liquid 
crystal hydrodynamics (||), (||), and (^ [Q. 



IV. PHASE ORDERING KINETICS 



A. Measures 

The kinetics of phase separation in liquid crystals has been examined using simulations by Zapotocky et al Q and 
more recently by Fukuda Q. The former studied phase separation in the diffusive regime whilst the latter added 
hydrodynamics. In both cases, rather deep quenches were performed and three different quantities considered in order 
to make a quantitative analysis of the domain growth. 

The first measure is the scalar correlation function for the tensorial nematic order parameter Q defined by |^ 

_ (Tr[Q(0,t)Q(r,t)]) 
(TrQ2(0,t)) ' 

where (• • •) denotes averaging over the positions 0. The correlation function is normalized so that C(0,t) = 1. A 
correlation length Lcor{t) at time t is defined by 
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C{L,or,t) = l/2. (21) 

Dynamical scaling states that the system is dynamically self-similar in time, except for a change in the length 
scale. If dynamical scaling holds, the correlation length will control the statistical properties of the system. Plotting 
the correlation function as a function of r/Lcor{t), the data at different times should collapse onto a single curve. 
Moreover Lcor{t) should decay with time as a power law 

L,or{t) ^ t^^"^ . (22) 

Zapotocky et al. Q obtained an exponent (j)cor — 0.41 significantly lower than the value 1/2 suggested by the diffusive 
character of the equation of motion for the order parameter (|[) and by scaling arguments ||2^ . 

A second measure is the Fourier transform of the correlation function, the structure factor S{k, t). The scaling form 
for the structure factor in d-dimensions is 

S{k,t)=Lt,{t)g{kL,or{t)) (23) 

where g should have the form giy) ~ y~(^+d) for large y for the 0{N) vector model (Porod's law) |2^. To see how 
this arises note that, for y = kLcor S> 1, the structure factor probes the order parameter at length scales much smaller 
than the separation between defects. Substantial variation of the order parameter on these length scales happens only 
in the vicinity of the defect cores, and is not related to inter-defect correlations. This implies that 

5(fc, t) - pdef{t)b{k), kL,or > 1, (24) 

where pdef is the density of defects in the system and h{k) is a function of k only. If we further assume that the 
separation of defects scales as the correlation length Lcor, so that pdef fx icor(0^^''^*'': where s is the dimensionality 
of the defect (zero here) then for d = 2 scaling implies 

9{V) - = y-\ (25) 

Zapotocky and collaborators did find the y~'^ behavior in the tails They also noted that the scaling functions 
appeared to approach the j/^^ law from above. This effect has been observed experimentally and it is attributed to 
inter-defect correlations. 

A scaling form also holds for the elastic energy. To see this note from Eq. (mj that 



J^el OC / dr{daQl3f )idaQl3j) 

dkk^S{k,t) 

= I dk[kL^or{t)]^g{kL^orit)) (26) 



for (i = 2 where the last line has been obtained by using the scaling form Eq. (|23|). The integral can be split into 
three parts, 0<r<a, a<r< Lcor{t) and Lcor{t) < r < cxj, where a is a cutoff. Assuming that the result is not 
dependent on the cutoff, the first part can be neglected, whereas for the second contribution using Porod's law for g 
and integrating from a to Lcor{t) gives 

Tel OC Leor{ty^ \n[LeoAt)/a]. (27) 

The last term can be neglected since it should scale roughly as L'^^^t) assuming that g{y) is approximately constant 
at small y. Differentiating Eq. ( |27| ) 

dlnTei d\nLcor[t) dXnTei < u\ ( n , ^ \ ^oo^ 

r{t)\-Z+—— , . , , (,2»j 



d In t dint d In Lcor \ \Ti[Lcor {t) I a 

where 0cor (i) is the exponent from the scaling of the correlation length. We define a new exponent for the characteristic 
energy as 

1 d In Tpj (f>cor (t) 
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Note the logarithmic corrections to the scaling of the elastic energy. Taking these into account, Zapotocky et al. 
obtained reasonable agreement between (jjgi and 0cor, whereas without taking them into account they obtained = 

0.325 <0eor §. 

A length scale can also be obtained from the average separation of topological defects. This separation is found by 
counting the number of defects in the system to obtain pdef and then taking 

Ldefit) = l/VAtef. (30) 

Zapotocky et al. obtained a value (j>def = 0.374 for the corresponding exponent. It is assumed that defect annihilation 
is the process which controls the growth of the correlation length and therefore dynamical scaling implies (jjcor = (f'def ■ 
Zapotocky et al. could give no explanation of the discrepancy in their values and interpreted it as an indication of 
the violation of dynamical scaling. 

Indeed in almost all of the published numerical work on phase ordering in systems which order via the annihilation 
of topological defects, it has been found that the exponents are less than the value 1/2 expected from dimensional 
analysis. A suggestive argument for the origin of this discrepancy has been given by Yurke et al § for the 0(2) model. 
These authors obtained an approximate equation of motion for an isolated defect-antidefect pair by equating the 
attractive and frictional forces acting on each defect. The attractive force was assumed to have the form Fat oc — 1/L 
where L is the separation of the defects. This would seem to be a reasonable assumption as it is well known the 



energy of such a pair goes like InL 17|. The frictional force was taken to be Ffr (x vln{R/a) where R is the "size" 



of a defect, a its core size, and = 5^ its velocity. The "size" is then assumed to be equal to L. Assigning a defect 
size of L may seem somewhat questionable. However the assumption can be supported by a more rigorous argument 

0- . . . 

Equating the frictional and elastic forces gives an implicit formula for the defect separation as a function of the 
time t before annihilation 



Ldef{t) 



t 



\n[Ldef{t)/a] - 1/2 



1/2 



(31) 



The growth in this length scale for increasing times before annihilation is then assumed to be the same as the growth 
in the average separation of defects after a quench. This gives an apparent exponent 



(l)def OC i 



t 



Llefit) HLdef{t)/a] 



(32) 



The argument implies that the failure to measure t^^^ is due to logarithmic corrections. However it has been claimed 
1^ that this is effectively untestable because of the unknown constant of proportionality in Eq. (^). One can however, 
measure an effective exponent for the two defect case and compare to that for the average separation after a quench. 
Zapotocky et al obtained a similar value, (j>def — 0.375 ^. 

More recently Fukuda has investigated the effect of the stress-induced flow on the kinetics of the nematic- 
isotropic transition by numerically solving the hydrodynamic equations for the tensor order parameter Q and the 
fluid velocity in c? = 2. The equations of motion were slightly different from the ones considered here. An older 
model which considers a corotational derivative rather than an upper convective derivative j2^] was used. For the 
purely dissipative case Fukuda obtained similar exponents to Zapotocky et al; for the hydrodynamic case he obtained 
(j)cor — 4'def = 0.43. Notc that, in contrast to the case without flow, 4>cor — 4'def indicating that the dynamical scaling 
hypothesis is confirmed when hydrodynamics is included. However the value calculated is still lower than the value 
1/2 expected for these systems. 

Most experiments on phase ordering in liquid crystals are in the three dimensional regime where the defect line 
energy controls the dynamics. There are only a limited number of experiments in a two dimensional geometry [^,^. 
Even in the experiments it is difficult to reach the asymptotic regime and the expected power law of for the 
number of defects is only obtained for the last decade of time m . 



B. Corrections to Scaling and Parameter space 



Our goal is to ascertain why previous studies have not measured the expected exponents and then, if possible, to 
perform simulations in the asymptotic scaling regime to obtain precise exponent values. In order to do this we first 
determine numerically the dependence of the corrections to scaling on the parameters at our disposal. In doing this 
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we will see that the difSculty in obtaining the asymptotic values of the exponents may have had as much to do with 
the fine scale details of the simulation on the lattice as with the constraints of lattice size and run times. 

As it seems clear from all previous studies that the phase behavior is dominated by the motion of Q in the plane, 
we restrict the order parameter to be of the form 

(Qxx Qxy \ 

Qxy Qyy , (33) 

~iQxx + Qyy) J 

in order to speed up the numerics. The quenches that we consider are envisaged to start at a high temperature, 
at which the equilibrium phase is the disordered, isotropic phase. We take as the initial condition for the lattice 
Boltzmann simulations a configuration representative of this phase in which the Qoif3 at each lattice point are random 
numbers uniformly distributed in a small interval [—6, S]. More precisely, by writing 

Qxx = A{3 cos (j) cos (j) — l)/2, 
Qxy = 3A cos sin 0/2, 

Qyy = A{3 sin (j) sin 2 (34) 

the initial configuration is obtained by assigning, at each lattice point, a random number between [—0.02,0.02] for 
the amplitude A and a random number between [— tt, tt] for the phase (j). 

To ascertain a rough dependence of the corrections to scaling on the various parameters of the system we examined 
quenches below and close to the spinodal line for several sets of parameters values. For this exploration of parameter 
space we used a lattice size of 256 x 256 and averaged over only two to five runs. 

Figure |l| shows the Schlieren pattern and associated fluid flow from a typical run. The shading indicates the 
orientation of the director field and the arrows the fluid velocity. The thing to note in such a diagram is the 
intersection of light and dark "brushes" which indicate the location of the ±1/2 disclinations. It is apparent that the 
main features of the flow are the vortices associated with the moving disclinations. The interaction of vortices in two 
dimensional flow has the same form as the interaction of disclinations in the director field so one does not expect any 
qualitative changes in the scaling behavior. 

We now compare 5 sets of data. The first four data sets are a quench to 7 = 3.25 (note that the first order transition 
is at 7 = 2.7 and that the spinodal is at 7 = 3.0). The first is a baseline run. In the second T is given by Eq.(|^) 
whereas in the others T is taken to be a constant and equal to 1 . In the third we turn off hydrodynamics by relaxing 
the momentum conservation constraint (i.e. the second Eq. in (p^)) at time t = 5000. The fourth has a value of k 
ten times smaller than the rest. The fifth data set is similar to the first except now the quench is deeper, to 7 = 3.7. 

The correlation function for 7 — 3.25 and with the renormalized diffusion constant is shown in Figure || There is 
a reasonable collapse of the data when the correlation function is plotted against r/Lcor{t)- This is as expected from 
dynamical scaling. However even if a function does appear to collapse it is very hard to measure the quality of the 
collapse and hence this provides only a weak test of scaling. 

Figure ^ shows the length scale obtained from the correlation function L^or a function of time for all five data 
sets. The exponents ^con 4'defj <t'ei obtained by fitting this data, as well as similar data for the number of defects and 
the distortion energy, to power laws are given in Table |. Although, as there are corrections to scaling, this is not a 
good way of obtaining the exponents, it does allow for a clear comparison between the different parameter sets. A 
more sophisticated way of estimating the growth exponents is to obtain an effective exponent (j){t*) by performing a 
simple linear fit to the plots of InL versus hii (e.g. Figure ||) for data in the interval i-m, i-m+i ■ ■ ■ ,t* ■ ■ tm-i,tm- 
TO = 8 is used in the fits presented here. 

From these preliminary runs with L = 256 we can make the following observations: 

• A renormalized diffusion constant F seems to help in the sense that the exponents are closer to the expected 
asymptotic value of 1/2. However this turns out to be mostly an early time effect. To see this we plot the effective 
correlation length exponent as a function of time in Figure ^ for the first two data sets in Table |. We see that for 
the baseline data set, the exponent first decreases and then increases with time whereas with the renormalized 
diffusion constant F the exponent immediately starts increasing with time toward the asymptotic value. Thus, 
we are effectively accessing longer times with the renormalized diffusion constant. This is reasonable because 
the initial ordering in the system, where the order parameter saturates away from the defect cores, is a diffusive 
process. 

• Hydrodynamics appears to speed the ordering very slightly, but not by a significant amount. The crucial test 
here is the run denoted in Figure ^ by filled triangles where at a given time step (t=5000) we have explicitly 
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switched off the hydrodynamics by randomizing the velocities directions after each colhsion step. One can see 
in Figure ^ that no significant change occurs with respect to the run in which the hydrodynamics is fuUy taken 
into account. This conclusion agrees with previous results by Fukuda even though his formulation of the 
model is somewhat different from ours. 

• Reducing k reduces the size of the defect cores and we find that changing k has a significant effect on the results. 
To see this we compare results for the baseline run with k = 0.02 to the run with a value of k ten times smaller 
(k — 0.002). The effective exponent extracted from the correlation function is shown in Figure ||. There are a 
few important things to note from this Figure. First the effective exponent is not a monotonic function of time: 
it reaches a maximum value at a given time and then it decreases as time is further increased. This is true for 
the baseline case as well, but the time at which it happens is much larger than for the run with more localized 
defects. Note that this effect has also been observed in previous simulations and interpreted as the onset of 
freezing in the dynamics, due mainly to the finite size of the system simulated. However the fact that the effect 
is stronger for smaller values of k suggests instead a lattice discretization effect. 

As an obvious length scale in this problem is the size of a defect core a it is tempting to take k as small as possible 
to minimize this scale. However, once a defect core is of the order of a lattice spacing, it becomes energetically 
advantageous for a defect to be centered in a plaquette of the lattice as this will minimize the distortion in the 
amplitude of the order parameter which is actually realized on the lattice, as shown in Figure |^. This creates an 
additional potential which can trap the defect. The height of this potential is a decreasing function of k as, if 
the amplitude distortion caused by a defect is spread over many sites, the defect will be less sensitive to exactly 
where it is centered. 

This effect occurs at later times in the simulation because inter-defect interactions decrease as the defect spacing 
increases. At the point where the lattice pinning potential balances the defect-defect interactions motion will 
stop. Indeed even before this time, the lattice pinning will cause an additional frictional drag on the defects 
which could be detrimental in trying to fit to theory. Thus we see that taking k small to try to achieve longer 
length scales (as appears to have been done in almost all previous work) is a trap better avoided. 

Second, as can be seen from Figure || , if the time is scaled by k the effective exponent for the baseline case 
appears to continue on the same curve as for smaller k. At first this seems surprising since one might expect 
time to scale with the diffusion constant, but not necessarily with n. To see why it occurs one can substitute 
a form for the director of a uniaxial nematic away from the defect core, n = (cos f?, sin 0) into Equation (|^). 
Ignoring the hydrodynamic flow one obtains a diffusion equation for 9 

dt9 = DeV^e, (35) 

where Dg = , confirming that the relevant diffusion constant is proportional to both f and k [ p^ , pl| . 

• Finally, the data for deeper quenches show an exponent further from 1/2. This is also due to the pinning 
potential described above. As can be seen from the effective exponents shown in Figure |^ the freezing happens 
at an earlier time for the deeper quench. This is not surprising as defects are more localized further from the 
transition and hence the lattice pinning potential will be greater. 



V. ASYMPTOTIC RESULTS: CRITICAL EXPONENTS 



In the previous section we examined the kinetics of phase ordering in a system of linear size L = 256 averaging over 
only a few initial configurations. This allowed us to explore in a more systematic way the space of the parameters 
(i.e. 7, K, r) relevant to the kinetics. In particular we found that by using a renormalized diffusion constant F and 
values of k big enough to avoid pinning effects better scaling behavior is achieved in less time. Using this knowledge 
we now simulate bigger systems and average over more starting configurations to obtain more precise results for the 
dynamical critical exponents of the model. 

We focus on the scaling behavior of the correlation function, the average defect separation and the elastic energy 
for phase separating systems of linear size L = 512. Averages were taken over 20 initial configurations. Quenches 
were run for 7 = 3.25 and 7 = 3.0. In both cases k = 0.3 and f = 1/(1 - ^TrQ'^)^. 

Figure ^ shows the time evolution of the length scale obtained from the correlation function for the two data sets. 
Note that at late times the data (on the log-log plot) displays a linear behavior with slope close to 1/2. A simple 
linear fit of the data with t > = 2239 gives 

0cor = 0.50 ±0.01 (36) 
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for 7 = 3.25, and 



(/)cor = 0.48 ± 0.02 (37) 

for 7 = 3.00. The errors were obtained by varying between 502 and 5012. This systematic error was much more 
significant than statistical errors evaluated from the statistics of the linear fit. 

As we have argued in the previous section a better quantity to check is the effective exponent as a function of 

fitting time t* . In Figure ^ we show the effective exponent for the correlation length (j)cor {t* ) as a function of inverse 
time. Note that the asymptotic value of 1/2 is clearly reached for the last decade in time. As far as we know this is the 
first time that strong evidence has been obtained for the conjectured 1/2 dynamical exponent for a phase separating 
nematic liquid crystal. 

Similar plots for the effective exponent for the number of defects (j)def{t*) are shown in Figure By looking at 
the high t* behavior we obtain 

(Ij^^j: = 0.47 ± 0.03 for 7 = 3.25, c/jdef = 0.48 ± 0.03 for 7 = 3.00 (38) 

where the errors have been estimated from the fluctuations in the effective exponents from a moving average with 
fewer points ( 4 in this case). Again, we see that the asymptotic value of the exponent is consistent with 1/2. 

Finally Figure |ll|(a) shows the exponent obtained from the elastic energy (j)ei{t*) . In this case a simple linear fit 
produces an asymptotic value smaller than 1/2 (~ 0.4 in the plot). This is mainly due to the logarithmic corrections 
to the scaling of the elastic energy (Eq. ^) that depress the effective growth exponent for the characteristic energy 



length with respect to the correlation length exponent (Eq. 29). Indeed by adding the logarithmic corrections to 
scaling as in Equation we can bring the exponent into agreement with that obtained from the correlation function, 
as shown in Figure pT|(b). 

By extrapolating the effective exponents at high t* values we obtain: 

(j)ei = 0.47 ± 0.03 for 7 = 3.25, 0ei = 0.50 ± 0.03 for 7 = 3.00. (39) 



VI. SHALLOW QUENCHES 

So far, we have examined the late stage kinetics of the orientational ordering, long after the amplitude has ordered. 
In this regime, the ordering has shown dynamics consistent with that expected for a diffusive XY model. It is also 
interesting to examine the early stage dynamics and shallower quenches where one may observe an interaction between 
the two symmetries of the problem (i.e. those of the scalar order parameter and the essentially 0(2) symmetry of the 
orientational order). 

First we examine what is naively expected for the amplitude ordering dynamics. Assuming that the orientational 
degrees of freedom are completely ordered the free energy can be written in terms of just the amplitude of the order 
parameter q 

+ (40) 

This is shown as the solid curves in Figure |lj for the cases 7 — 2.7 (phase transition point) and 7 — 3.0 (spinodal 
line, where the barrier between the isotropic and nematic states vanishes). One expects the presence of a free energy 
barrier between the isotropic and nematic states to affect the phase ordering dynamics. In the presence of a barrier 
a finite region of the nematic phase must be nucleated to initiate growth. Beyond the spinodal (which occurs at 
7 = 3.0) the system will order spontaneously by spinodal decomposition throughout space. If we artificially orient 
the director field before quenching, this is indeed what is observed. 

However, in a normal quench where the director is not oriented, one observes somewhat different behavior. In 
particular, there is a region where one might expect spontaneous ordering (3.0 < 7 < 7s (k, L)) which in fact exhibits 
a two stage growth process in the absence of a nucleating site, and domain growth from a nucleating site (as opposed 
to spinodal decomposition) if one is provided. 

Figure |l^ shows the amplitude as a function of time for 7 = 3.05 for a system which is initially disordered and 
in which we have not supplied any strong nucleating sites. There are clearly two different stages of growth. In the 
first stage the amplitude grows very slowly for a long time. Then the system abruptly orders, very quickly reaching 
the equilibrium value. If instead we supply the system with nucleating sites at the beginning of the simulation it 
orders much more quickly with domains growing from the nucleating sites. These observations are consistent with 
the presence of a free energy barrier which decays with time. 
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The free energy barrier can be identified as resulting from the free energy of the orientational order. There is a 
time dependent quadratic term in the free energy due to the distortion. The distortion energy is proportional to nq^ 
(see Eq. (|l|)) so that one has an effective term 

Td cx Kq^Ldefit)-^ \ll[Ldef{t)/a], (41) 



using Eq. (27) with Ldef{t) being the separation between defects. As a result, not only does the state at 7 = 3.0 
remain metastable but the global free energy minimum at q > may not be a minimum at zero time. This contribution 
to the free energy is included in Figure n3 for different values of Ldef ■ 

The domain growth in the nucleated case is now easily understood. Why there is eventually ordering with no nucle- 
ation sites still requires some explanation. Eq. (|4^) gives the average free energy contribution from the orientational 
order. However, locally there are regions where there is higher orientational order so that some gain in free energy can 
be made by increasing q in these regions. These regions are still not large enough to gain sufRciently from this to pay 
for the cost of an interface between fully ordered and disordered states so the elastic energy suppresses the growth in 
q, keeping its typical value very small. However, as we have already pointed out in Eq. (^5|) the diffusion constant 
for the orientational ordering is independent of q, as long as g > 0. Hence the orientational ordering can proceed. 
(See, for instance, the number of defects as a hmction of time shown in Figure |lj.) As the orientation proceeds, Ldef 
grows and eventually the free energy barrier is small enough that the system can form a sufficiently large region of 
the amplitude ordered phase to nucleate growth. Once such a region is formed the amplitude ordering can proceed 
rapidly. 



VII. SUMMARY AND DISCUSSION 



To summarize: we have used a lattice Boltzmann algorithm to investigate the phase ordering of liquid crystals 
quenched from the isotropic to the nematic phase. The system is described by the Beris-Edwards equations of liquid 
crystal hydrodynamics. These are written in terms of a tensor order parameter which means that the dynamics of 
topological defects which drives the phase ordering appears naturally in the simulations. The liquid crystal moves 
towards an equilibrium described by the Landau-de Gennes free energy. 

The flexibility of the numerical scheme allows an investigation into how the ordering process is affected by the 
model parameters. In particular we find that if the length scale defined by the size of the topological defects is taken 
to be too short the cores are pinned by the lattice as the simulation proceeds giving an incorrect value for the growth 
exponent. This is likely to have been the problem faced by other authors who consistently obtained a value less than 
the 1/2 expected from scaling arguments. 

Averaging over 20 runs on a lattice of linear size 512 we obtained 0cor = 0.49 ± 0.02, = 0.485 ± 0.03, (pdef — 
0.475±0.03. The value for was fit using logarithmic corrections. Hydrodynamic flow speeds up the ordering slightly 
but does not change the scaling behavior. This is reasonable because, in two dimensions, interactions between the 
flow vortices set up by the moving defects have the same logarithmic form as the defect-defect interactions themselves. 
Thus theoretical and numerical results arc now in pleasing agreement. 
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run 


parameter set 




exponents 

<t>def 


0d 


baseline 


7 = 3.25, r = 1,K = 0.02 


0.42 


0.406 


0.35 


renorm. diffusion 


7 = 3.25, r = (1 - |rrQ2)-^ k = 0.02 


0.45 


0.43 


0.365 


no hydrodynamics 


7 = 3.25, r = Ik = 0.02 


0.42 


0.398 


0.345 


more localized defects 


7 = 3.25, r = 1,K = 0.002 


0.36 


0.33 


0.27 


deeper quench 


7 = 3.7, r = 1,K = 0.02 


0.39 


0.36 


0.31 



TABLE I. Comparison of the exponents obtained by fitting the data from runs on small systems (e.g. Fig. ^ to a power law. 
The statistical errors ^ 0.003 and the systematic errors ^ 0.015 (as the end-points of the fit are changed). 
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FIG. 1. Schlicrcn pattern (shading) and fluid velocity (vectors) associated with the ordering of a liquid crystal after a quench 
from the isotropic to the nematic phase. The director in the darkest regions is perpendicular to the director in the lightest 
regions (with the other shade in between). The shading on the Schlieren pattern has been rendered with only three shades of 
gray so as not to obscure the vector field. Disclinations of strength ±1/2 are located at the intersection of dark grey and white 
"brushes". 
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FIG. 2. Correlation function for 7 = 3.25 and with renormalized diffusion constant. Inset: unsealed correlation function for 
times between t = 159 and 8913. Main figure: correlation function as a function of r/L{t) where L{t) is determined as the 
length at which the unsealed correlation function crosses 1/2. 
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FIG. 3. Length scale obtained from sealing of correlation function as a function of time. The lines are power law fits and 
the symbols correspond to box: 7 = 3.25 baseline, diamond: 7 = 3.25 renormalized diffusion constant, filled triangle: 7 = 3.25 
hydrodynamics off, star: 7 = 3.25 small kappa, hollow triangle: 7 = 3.7 as baseline. In this figure (and all others) times are 
multiplied by F to make them dimensionless and distances are measured in units of the lattice spacing of the simulation. 
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FIG. 4. Effective exponent for the growth of the correlation length as a function of time empty triangle: 7 = 3.25 baseline, 
filled triangle : 7 = 3.25 renormalized diffusion constant. 
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FIG. 5. Effective exponent for the growth of the correlation length as a function of time empty triangle: 7 = 3.25 baseline, 
filled triangle: 7 = 3.25 k = 1/10 of the baseline value. The time t* has been multiplied by 1/10 (i.e. time is scaled by 1/k) 
for the small kappa case. 
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FIG. 6. Two largest eigenvalues of Q along a cut joining the center of two ±1/2 defects moving to annihilate each other. 
The squares and diamonds are for the case with k — 0.02 and the sta rs and triangles for n = 0.5. The dips are an indication 
of the lattice pinning potential. The abscissa is scaled by the length -y/ k/pq to make it dimensionless (all other parameters are 
the same). 
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FIG. 7. Effective exponent for the growth of the correlation length as a function of time triangle: 7 = 3.25 baseline, box: 
7 = 3.7 baseline. 
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FIG. 8. Length scale obtained from the correlation function as a function of time. The eyeguide line has slope 1/2. The 
symbols correspond to empty triangle: 7 = 3.25 and filled triangle: 7 = 3.00. 
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FIG. 9. Effective exponent for the growth of the correlation length as a function of time: filled triangle: 7 = 3.00, empty 
triangle: 7 = 3.25. 
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FIG. 10. Effective exponent for the growth of the defect separation length as a function of time filled triangle: 7 = 3.00, 
empty triangle: 7 = 3.25. 
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FIG. 11. (a)Effective exponent for the growth of the elastic energy length as a function of time triangle: 7 — 3.00, box: 
7 — 3.25. (b) Effective exponent for the correlation length (filled triangle), and the growth of the elastic energy corrected for 
the logarithmic scaling of Eq.(^9|) (empty triangle). 
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FIG. 12. Free energy (xlO"') versus magnitude of the order parameter q for 7 = 2.7 (left) and 7 = 3.0 (right). Solid lines 
correspond to the assumption of complete orientational order and dashed lines correspond to the spatially averaged free energy 
for different values of the defect separation (6 and 12 in left figure and 2 in right figure). 
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FIG. 13. Amplitude q versus time for a quench to 7 = 3.05. The initial fractional noise on q (~ 10 ®) is very small and 
cannot cause ordering by nucleation and growth. 
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